# M. Orsi - Essex group - October 2006
# modified December 2006, to consider total area + volume
# modified december 2008, to deal with only two columns (area + volume) - got rid of first column (time) to save on disk space
# This program opens the file "areaVol.dat", which is one of BRAHMS's output files, and extract the area compressibility (AKA stretch) modulus

# USAGE: cat area_volume.dat | python computeAreaCompressibilityModulus.py
         
import sys, string
from math import sqrt
lines = sys.stdin.readlines()

#constants
BOLTZMANN = 1.3806505e-23 # [J / K]

# variables
T = 303.15 # [K]
nLipids = 128 # number of lipid molecules in the system

print "REMEMBER TO ADJUST THE NUMBER OF LIPID AND TEMPERATURE!!!"

lineCounter = 0
meanLipidAreaSum = 0
squareLipidAreaSum = 0
meanBoxAreaSum = 0
squareBoxAreaSum = 0

brahmsAreaVolFile = open( 'area_volume.dat', 'r' )

for line in lines: # only consider boxArea!
    lineCounter = lineCounter + 1
    words = string.split( line )
    boxArea = string.atof(words[0]);
    meanBoxAreaSum = meanBoxAreaSum + boxArea;
    squareBoxAreaSum = squareBoxAreaSum + boxArea**2

brahmsAreaVolFile.close

# computing statistical quantities

meanBoxArea = meanBoxAreaSum / lineCounter;
print "meanBoxArea = %4.2f A^2" % (meanBoxArea * 100)
meanSquaredBoxAreaFluctuation = squareBoxAreaSum / lineCounter - meanBoxArea**2
# print "meanSquaredAreaFluctuation = %6.3f A^4" % (meanSquaredAreaFluctuation * 10000)

meanLipidArea = meanBoxArea / (nLipids/2);
print "meanLipidArea = %4.2f A^2" % (meanLipidArea * 100)

# computing modulus
KA_boxArea = BOLTZMANN * T * meanBoxArea / meanSquaredBoxAreaFluctuation # [ J / nm^2 ]

# conversion considering that J/nm^2 = 10^21 dyn/cm
KA_boxArea = KA_boxArea * pow(10,21)

print "KA = %6.1f dyn/cm" % ( KA_boxArea )
